Figure 1: code for insets
Coverage example
## Project population
vacc_reconstruct <- function(t, y, parms){
dS <- -parms["deaths"] * y["S"] + parms["births"] * (y["S"] + y["V"]) + parms["waning"] * y["V"]
dV <- -parms["deaths"] * y["V"] - parms["waning"] * y["V"]
dWaning <- -parms["waning"] * y["V"]
dDeaths <- -parms["deaths"] * y["V"]
dBirths <- -parms["births"] * (y["S"] + y["V"])
return(list(c(dS, dV)))
}
eventfun <- function(t, y, parms){
with (as.list(y),{
diff <- 0.7*(S + V) - V ## newly vaccinated
V <- V + diff
S <- S - diff
return(c(S, V))
})
}
## times in relevant timestep
times <- seq (1/52, 5, by = 1/52) ## weekly for 10 years
## initial pop
IS <- 70000
IV <- 0
## param
births <- 0.50
deaths <- 0.42
waning <- 0.33
parms <- c(deaths = deaths, births = births, waning = waning)
## run
calcsus <- lsoda(y = c(S = IS, V = IV), times = times, func = vacc_reconstruct, parms = parms,
events = list(func = eventfun, time = seq(1/52, 5, by = 1)))
vacc_df <- data.frame(Year = calcsus[, "time"], propVacc = calcsus[,"V"]/(calcsus[,"V"] + calcsus[,"S"]))
vacc_a <- ggplot(data = vacc_df, aes(x = Year, y = propVacc)) +
geom_line(color = "lightseagreen", size = 1) +
xlab ("Year") +
ylab("Coverage") +
ylim(c(0, 1)) +
geom_hline(yintercept = 0.69, color = "aquamarine4", linetype = 2, size = 1) +
newtheme +
theme(text = element_text(color = "aquamarine4"), axis.text = element_text(color = "aquamarine4"),
axis.line = element_line(color = "aquamarine4"))
ggsave("figs/fig1/vacc_a.jpeg", vacc_a, device = "jpeg", width = 5, height = 5)
# include_graphics("figs/fig1/vacc_a.jpeg")
Secondary cases
## params
secondaries <- data.frame(secondaries = rnbinom(1000, mu = 1.2, size = 1.3))
sec_b <- ggplot(data = secondaries, aes(x = secondaries)) +
geom_histogram(binwidth = 1, color = "lightgray", fill = "red", alpha = 0.5) +
geom_vline(xintercept = 1.2, color = "red", linetype = 2, size = 1.1) +
xlab("Secondary cases") +
ylab("Frequency") +
newtheme +
theme(text = element_text(color = alpha("firebrick3", 1)),
axis.text = element_text(color = alpha("firebrick3", 1)),
axis.line = element_line(color = alpha("firebrick3", 1)))
ggsave("figs/fig1/sec_b.jpeg", sec_b, device = "jpeg", width = 5, height = 5)
# include_graphics("figs/fig1/secc_b.jpeg")
Incubation period
incubation <- as.data.frame(list(Days = seq(0, 365, 1),
Density = dgamma(seq(0, 365, 1), shape = 1.1518,
rate = 0.0429)))
inc_c <- ggplot(data = incubation, aes (x = Days, y = Density)) +
geom_line(color = "darkred", size = 1.2) +
xlab("Incubation Period \n (days)") +
geom_vline(xintercept = 22.3, color = "darkred", linetype = 2, size = 1.1, alpha = 0.5) +
newtheme +
theme(text = element_text(color = "darkred"), axis.text = element_text(color = "darkred"),
axis.line = element_line(color = "darkred"))
ggsave("figs/fig1/inc_c.jpeg", inc_c, device = "jpeg", width = 5, height = 5)
# include_graphics("figs/fig1/inc_c.jpeg")
Infectious period
infectious <- as.data.frame(list(Days = seq(0, 10, 0.1),
Density = dgamma(seq(0, 10, 0.1), shape = 2.9012,
rate = 1.013)))
inf_d <- ggplot(data = infectious, aes (x = Days, y = Density)) +
geom_line(color = "navy", size = 1.2) +
xlab("Infectious Period \n (days)") +
geom_vline(xintercept = 3.1, color = "navy", linetype = 2, size = 1.2, alpha = 0.5) +
newtheme +
theme(text = element_text(color = "navy"), axis.text = element_text(color = "navy"),
axis.line = element_line(color = "navy"))
ggsave("figs/fig1/inf_d.jpeg", inf_d, device = "jpeg", width = 5, height = 5)
# include_graphics("figs/fig1/inf_d.jpeg")
Road network with scale
tz_roads <- readOGR("data/Tanzania_Roads/Tanzania_Roads.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/mrajeev/Documents/Projects/ModelingChapter/data/Tanzania_Roads/Tanzania_Roads.shp", layer: "Tanzania_Roads"
## with 518 features
## It has 17 fields
tz_roads <- fortify(tz_roads)
map_e1 <- ggplot(data = tz_roads, aes(x = long, y = lat, group = group)) +
geom_line(color = "mediumorchid4") +
scalebar(tz_roads, location = "bottomleft",
dist = 200, dist_unit = "km", transform = TRUE, box.fill = c("mediumorchid4", "grey"),
box.color = c("mediumorchid4", "grey"),
st.size = 5, st.dist = 0.05, st.color = "mediumorchid4",
height = 0.02, model = 'WGS84') +
theme_map()
ggsave("figs/fig1/map_e1.jpeg", map_e1, device = "jpeg", width = 5, height = 5)
# include_graphics("figs/fig1/map_e1.jpeg")
Dispersal kernel
## params
dispersal <- as.data.frame(list(km = seq(0, 10, 0.01),
Density = dgamma(seq(0, 10, 0.01), shape = 0.215,
scale = 0.245)))
disp_e2 <- ggplot(data = dispersal, aes (x = km, y = Density)) +
geom_line(color = "mediumorchid4", size = 1.2) +
xlab("Distance to next bite \n (km)") +
geom_vline(xintercept = 0.88, color = "mediumorchid4", linetype = 2, size = 1.2, alpha = 0.75) +
newtheme +
theme(text = element_text(color = "mediumorchid4"), axis.text = element_text(color = "mediumorchid4"),
axis.line = element_line(color = "mediumorchid4"))
ggsave("figs/fig1/disp_e2.jpeg", disp_e2, device = "jpeg", width = 5, height = 5)
# include_graphics("figs/fig1/disp_e2.jpeg")
Figure 2: Summarizing studies
studies <- read.csv("data/modeling_studies.csv")
## A. Studies by Year
year_A <- ggplot(data = studies, aes(x = Year)) +
geom_bar(fill = "grey50") +
scale_y_continuous(labels = scales::number_format(accuracy = 1)) +
ylab("Number of studies") +
labs(tag = "A") +
newtheme +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
## B. Studies by country
## refactor countries first
studies %>%
mutate(Country = fct_recode(Country, Other = "India", Other = "Israel", Other = "Kenya")) -> studies
## then plot
country_B <- ggplot(data = studies, aes(x = fct_relevel(reorder(Country, Country, function(x)-length(x)),
"Multiple (Africa)", "Multiple (global)",
"Other", "Hypothetical", after = 10))) +
geom_bar(fill = "grey50") +
scale_y_continuous(breaks = seq(0, 10, by = 2), labels = scales::number_format(accuracy = 1)) +
ylab("Number of studies") +
xlab("") +
labs(tag = "B") +
newtheme +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
## C. R0 plot
r_ests <- read.csv("data/r_ests.csv")
R0_C <- ggplot(data = filter(r_ests, R_type == "R0"), aes(x = Estimate)) +
geom_histogram(binwidth = 0.25, color = "white", fill = "grey50") +
scale_y_continuous(labels = scales::number_format(accuracy = 1)) +
labs(tag = "C") +
xlab("R0 estimate") +
ylab("Number of studies") +
newtheme
## D. Key modeling decisions and features
studies %>%
select(14:20) %>%
gather() %>%
filter(value %in% "Y") %>%
group_by(key) %>%
summarize(proportion = n()/51) -> summary
features_D <- ggplot(data = summary, aes(x = reorder(key, -proportion), y = proportion)) +
geom_col(fill = "grey50") +
scale_x_discrete(labels = c("DDT",
"Fit to data", "Stochastic", "Spatial",
"Heterogeneity", "Introductions",
"Underreporting")) +
ylab("Proportion of studies") +
ylim(c(0, 0.8)) +
xlab("") +
labs(tag = "D") +
newtheme +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
fig2 <- year_A + R0_C + country_B + features_D + plot_layout(ncol = 2, nrow = 2)
ggsave("figs/fig2.jpeg", fig2, device = "jpeg", width = 20, height = 15)
## Warning: Removed 1 rows containing non-finite values (stat_count).
include_graphics("figs/fig2.jpeg")
Figure 3: data
## Reported data
repdata <- read.csv("data/reported_data.csv")
repdata %>%
mutate(Nobs_maxed = case_when(Nobs < 250 ~ as.numeric(Nobs), Nobs >= 250 ~ 250),
Spatial.scale = fct_relevel(Spatial.scale, "National", after = Inf),
Temporal.scale = fct_relevel(Temporal.scale, "Month", after = 2)) -> repdata
data_A <- ggplot(data = repdata, aes(x = reorder(Type.of.data, Type.of.data,function(x)-length(x)))) +
geom_bar() +
newtheme +
scale_x_discrete(labels = c("Clinical and lab \n confirmed animal cases",
"Reported human \n deaths",
"Sequence data", "Lab confirmed \n animal cases",
"Lab confirmed \n human exposures", "Reported \n animal bites")) +
xlab("Type of data") +
ylab("Number of data \n sources") +
labs(tag = "A") +
coord_flip()
data_B <- ggplot(data = repdata, aes(x = Temporal.scale, y = Length, fill = Spatial.scale, size = Nobs_maxed)) +
ggbeeswarm::geom_beeswarm(cex = 3.5, alpha = 0.75, shape = 21, color = "grey50", stroke = 1.5) +
scale_fill_manual(values = c("Darkred", "Red", "Pink", "Grey50"), name = "Spatial scale") +
scale_size_area(name = "# of Observations", breaks = c(10, 100, 150, 250),
labels = c("10", "100", "150", "250+")) +
ylab("Length of time \n series (years)") +
xlab("Temporal scale") +
newtheme +
labs(tag = "B") +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
guides(fill = guide_legend(override.aes = list(size = 3)))
fig3 <- data_A + data_B + plot_layout(ncol = 2, widths = c(1, 1))
ggsave("figs/fig3.jpeg", plot = fig3, device = "jpeg", height = 6, width = 15)
include_graphics("figs/fig3.jpeg")
Supplementary Figure S1
# models
SEI.rabies <- function(nu = 0.50, mu = 0.42, R0 = 1.2, sigma = (1/22.3*365),
gamma = (1/3.1*365), K = 20, START.S = 50000,
START.E = 0, START.I = 2,
START.S.density = 15, START.I.density = 0.01,
START.E.density = 0.00, years = 20, steps = 1/52, model = "frequency",...){
### Deterministic skeleton of models for rabies
require(deSolve)
times <- seq(from = 0, to = years, by = steps)
## Simple FDT model with complete disease-induced mortality
if (model == "frequency"){
START.N <- START.S + START.E + START.I
beta = (R0*gamma*(sigma + mu))/sigma
# This function models a time step for the SIR:
dx.dt.SIR <- function(t, y, parms) {
N <- y["S"] + y["E"] + y["I"]
S <- y["S"]
E <- y["E"]
I <- y["I"]
# Calculate the change in Susceptible
dS <- parms["nu"]*(S + E) - parms["mu"]*S - parms["beta"]*S*I/N
#Calculate the change in Exposed
dE <- parms["beta"]*S*I/N - parms["mu"]*E - parms["sigma"]*E
# Calculate the change in Infected
dI <- parms["sigma"]*E - parms["gamma"]*I
# Return a list with the changes in S, E, I, R at the current time step
return(list(c(dS, dE, dI)))
}
# Create the parameter vector
parms <- c(nu = nu, mu = mu, beta = beta, sigma = sigma, gamma = gamma)
inits <- c(S = START.S, E = START.E, I = START.I)
}
## Classical fox rabies model (widely used for dog rabies as well, adapted from Anderson et al. 1981)
if (model == "density"){
START.N.density <- START.S.density + START.I.density + START.E.density
# the old school rabies model
# This function models a time step for the SIR:
dx.dt.SIR <- function(t, y, parms) {
N <- y["S"] + y["I"] + y["E"]
S <- y["S"]
E <- y["E"]
I <- y["I"]
## Calculate dmort and beta here from parms and plug them in
dmort <- (parms["nu"] - parms["mu"])/parms["K"]
beta <- (parms["R0"]*(parms["sigma"] + parms["nu"])*(parms["gamma"] + parms ["nu"]))/(parms["sigma"]*parms["K"])
# Calculate the change in Susceptible
dS <- parms["nu"]*N - parms["mu"]*S - dmort*S - beta*S*I
#Calculate the change in Exposed
dE <- beta*S*I - parms["mu"]*E - dmort*N*E - parms["sigma"]*E
# Calculate the change in Infected
dI <- parms["sigma"]*E - parms["mu"]*I - dmort*N*I - parms["gamma"]*I
# Return a list with the changes in S, E, I, R at the current time step
return(list(c(dS, dE, dI)))
}
# Create the parameter vector
parms <- c(nu = nu, mu = mu, R0 = R0, K = K, sigma = sigma, gamma = gamma)
inits <- c(S = START.S.density, E = START.E.density, I = START.I.density)
}
# Run the ODE solver
SIR.output <- lsoda(y = inits,
times = times,
func = dx.dt.SIR,
parms = parms)
SIR.output <- as.data.frame(SIR.output)
SIR.output$N <- SIR.output$S + SIR.output$E + SIR.output$I
return (SIR.output)
}
## Density dependent
ddt_1.01 <- SEI.rabies(R0 = 1.01, years = 20, model = "density")
ddt_1.01 <- data.frame(infected = unname(tapply(ddt_1.01$I, (seq_along(ddt_1.01$I)-1) %/% 4, sum)[1:260]),
pop = ddt_1.01$N[seq(1, 20*52, 4)], R0 = 1.01, trans = "Density")
ddt_1.1 <- SEI.rabies(R0 = 1.1, years = 20, model = "density")
ddt_1.1 <- data.frame(infected = unname(tapply(ddt_1.1$I, (seq_along(ddt_1.1$I)-1) %/% 4, sum)[1:260]),
pop = ddt_1.1$N[seq(1, 20*52, 4)], R0 = 1.1, trans = "Density")
ddt_1.05 <- SEI.rabies(R0 = 1.05, years = 20, model = "density")
ddt_1.05 <- data.frame(infected = unname(tapply(ddt_1.05$I, (seq_along(ddt_1.05$I)-1) %/% 4, sum)[1:260]),
pop = ddt_1.05$N[seq(1, 20*52, 4)], R0 = 1.05, trans = "Density")
ddt <- bind_rows(ddt_1.01, ddt_1.1, ddt_1.05)
ddt$time <- 1:260
ddt_plot <- ggplot(data = ddt, aes(x = time, y = infected/pop, group = R0, color = as.factor(R0))) +
geom_line(size = 1) +
scale_color_manual(values = c("lightcoral", "red", "darkred"), name = "R0") +
xlab("Months") +
ylab("Incidence (monthly proportion \n of population infected)") +
newtheme +
labs(tag = "B", subtitle = "Density")
## Frequency dependent
fdt_1.01 <- SEI.rabies(R0 = 1.01, years = 20, model = "frequency")
fdt_1.01 <- data.frame(infected = unname(tapply(fdt_1.01$I, (seq_along(fdt_1.01$I)-1) %/% 4, sum)[1:260]),
pop = fdt_1.01$N[seq(1, 20*52, 4)], R0 = 1.01, trans = "Frequency")
fdt_1.1 <- SEI.rabies(R0 = 1.1, years = 20, model = "frequency")
fdt_1.1 <- data.frame(infected = unname(tapply(fdt_1.1$I, (seq_along(fdt_1.1$I)-1) %/% 4, sum)[1:260]),
pop = fdt_1.1$N[seq(1, 20*52, 4)], R0 = 1.1, trans = "Frequency")
fdt_1.05 <- SEI.rabies(R0 = 1.05, years = 20, model = "frequency")
fdt_1.05 <- data.frame(infected = unname(tapply(fdt_1.05$I, (seq_along(fdt_1.05$I)-1) %/% 4, sum)[1:260]),
pop = fdt_1.05$N[seq(1, 20*52, 4)], R0 = 1.05, trans = "Frequency")
fdt <- bind_rows(fdt_1.01, fdt_1.1, fdt_1.05)
fdt$time <- 1:260
fdt_plot <- ggplot(data = fdt, aes(x = time, y = infected/pop, group = R0, color = as.factor(R0))) +
geom_line(size = 1) +
scale_color_manual(values = c("lightcoral", "red", "darkred"), name = "R0") +
xlab("Months") +
ylab("Incidence (monthly proportion \n of population infected)") +
guides(colour = "none") +
newtheme +
labs(tag = "A", subtitle = "Frequency")
figS1 <- fdt_plot / ddt_plot
ggsave("figs/figS1.jpeg", plot = figS1, device = "jpeg", height = 11, width = 11)
include_graphics("figs/figS1.jpeg")
Supplementary Figure S2
## Study category
studies %>%
mutate_at(vars(starts_with("Category")),
function(x) ifelse(x == "N" | is.na(x), 0, 1)) %>%
select(Reference, starts_with("Category")) -> categories
names(categories)
## [1] "Reference"
## [2] "Category..explain.observed.patterns"
## [3] "Category..estimate.key.parameters"
## [4] "Category..study.control.measures"
## [5] "Category..identify.drivers.of.dynamics"
## [6] "Category..predict.future.trends"
names(categories) <- c("Reference", "Explain patterns",
"Estimate parameters", "Study control",
"Identify drivers", "Predict trends")
## Upset plot
jpeg("figs/figS2.jpeg", width = 500, height = 500)
upset(categories, nsets = 5, nintersects = NA,
sets.x.label = "Type of study", mainbar.y.label = "Number of studies", text.scale = 2)
dev.off()
## quartz_off_screen
## 2
## Bar plot
studies %>%
select(Reference, starts_with("Category")) %>%
gather() %>%
filter(value %in% "Y") -> cat_summary
## Warning: attributes are not identical across measure variables;
## they will be dropped
table(cat_summary$key)
##
## Category..estimate.key.parameters
## 15
## Category..explain.observed.patterns
## 12
## Category..identify.drivers.of.dynamics
## 10
## Category..predict.future.trends
## 1
## Category..study.control.measures
## 26
include_graphics("figs/figS2.jpeg")
sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] UpSetR_1.4.0 patchwork_1.0.0.9000 forcats_0.4.0
## [4] stringr_1.4.0 dplyr_0.8.3 purrr_0.3.2
## [7] readr_1.3.1 tidyr_1.0.0 tibble_2.1.3
## [10] tidyverse_1.2.1 ggsn_0.5.0 ggthemes_4.2.0
## [13] rgdal_1.4-6 sp_1.3-1 deSolve_1.24
## [16] ggplot2_3.2.1 knitr_1.25
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.1 jsonlite_1.6 modelr_0.1.5
## [4] assertthat_0.2.1 vipor_0.4.5 cellranger_1.1.0
## [7] yaml_2.2.0 pillar_1.4.2 backports_1.1.5
## [10] lattice_0.20-38 glue_1.3.1 digest_0.6.22
## [13] rvest_0.3.4 colorspace_1.4-1 htmltools_0.4.0
## [16] plyr_1.8.4 pkgconfig_2.0.3 broom_0.5.2
## [19] haven_2.1.1 scales_1.0.0 jpeg_0.1-8
## [22] ggmap_3.0.0 generics_0.0.2 ellipsis_0.3.0
## [25] withr_2.1.2 lazyeval_0.2.2 cli_1.1.0
## [28] magrittr_1.5 crayon_1.3.4 readxl_1.3.1
## [31] maptools_0.9-5 evaluate_0.14 nlme_3.1-141
## [34] xml2_1.2.2 foreign_0.8-72 class_7.3-15
## [37] beeswarm_0.2.3 tools_3.6.1 hms_0.5.1
## [40] RgoogleMaps_1.4.4 lifecycle_0.1.0 munsell_0.5.0
## [43] compiler_3.6.1 e1071_1.7-2 rlang_0.4.1
## [46] classInt_0.4-1 units_0.6-4 rstudioapi_0.10
## [49] rjson_0.2.20 bitops_1.0-6 labeling_0.3
## [52] rmarkdown_1.16 gtable_0.3.0 DBI_1.0.0
## [55] R6_2.4.0 gridExtra_2.3 lubridate_1.7.4
## [58] zeallot_0.1.0 KernSmooth_2.23-15 stringi_1.4.3
## [61] ggbeeswarm_0.6.0 Rcpp_1.0.2 vctrs_0.2.0
## [64] sf_0.8-0 png_0.1-7 tidyselect_0.2.5
## [67] xfun_0.10
